#include <iostream>
#include <cmath>
#include <vector.h>

using namespace std;

vector<double> f(vector<double> u, double mu)
{
    double u1 = u[0];
    vector<double> f(6);
    f[0] = u[3];
    f[1] = u[4];
    f[2] = u[5];
    f[3] = 2 * u[4] + u[0] - mu * (u[0] + mu - 1)/(pow((u[1]*u[1] + u[2]*u[2]+(u[0]+mu-1)*(u[0]+mu - 1)),1.5))-(1-mu)*(u[0] + mu)/(pow((u[1]*u[1]+u[2]*u[2]+(u[0]+mu)*(u[0]+mu)),1.5));
    f[4] = -2 * u[3] +u[1]-mu * u[1]/(pow((u[1]*u[1]+u[2]*u[2]+(u[0]+mu-1)*(u[0]+mu-1)),1.5))-(1-mu)*u[1]/(pow((u[1]*u[1]+u[2]*u[2]+(u[0]+mu)*(u[0]+mu)),1.5));
    f[5] = -mu * u[2]/(pow((u[1]*u[1]+u[2]*u[2]+(u[0]+mu-1)*(u[0]+mu-1)),1.5))-(1-miu)*u[2]/(pow((u[1]*u[1]+u[2]*u[2]+(u[0]+mu)*(u[0]+mu)),1.5));
    return f;			   

};

class TimeIntegrator
{
public:
    vritual vector<vector<double>> solve(double k, double mu, vector<double> u0, int step, int p) = 0;
};

class ABM : public TimeIntegrator
{
public:
    vector<vector<double>> solve(double k, double mu, vector<double> u0, int step, int p)
	{
	    vector<vector<double>> U;
	    if (p == 1)
	    {
		vector<double> u(6), t(6);
		u = u0;
		U.push_back(u0);
		for (int i = 1; i <= step; i++)
		{
		    u = u + k * f(u, mu);
		    t = u;
		    U.push_back(t);
		}
	    }else if (p == 2)
	    {
		vector<double> u(6), u_1(6), t(6);
		u_1 = u0;
		U.push_back(u0);
		u = u_1 + k*f(u_1, mu);
		t = u;
		U.push_back(t);
		for (int i = 1; i <= step - 1; i++)
		{
		    vector<double> tmp(6);
		    tmp = u;
		    u = u + k*(1.5*f(u, mu) - 0.5*f(u_1, mu));
		    u_1 = tmp;
		    t = u;
		    U.push_back(t);
		}
	    }
	};
};
